//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

#define WIDTH 800
#define HEIGHT 800

#define COILMATRIX_SIDE 13

int main (int argc, char *argv[])
{
  prec_t radius = 0.1524;
  prec_t coil_side = 0.0508;
  prec_t thickness = prec_t_sqrt (2 * coil_side * coil_side);
  prec_t wr = thickness / 2;
  prec_t spacing = 0.01;
  prec_t rr = radius + (wr * 2 + spacing) / prec_t_sqrt(2.0);
  prec_t I = 134560;

  vect3d c1(-rr, radius, 0);
  vect3d c2(-radius, rr, 0);
  vect3d c3;
  vect3d mid = 0.5 * (c1 + c2);

  prec_t fmin, fmax;

  ShapedCoilFactory round_shaper(COILMATRIX_SIDE);
  round_shaper.set_circle();
  // divide by sqrt2 so would fit withing circle of wr
  ShapedCoilFactory rect_shaper(COILMATRIX_SIDE, COILMATRIX_SIDE, 1/PREC_SQRT2, 1/PREC_SQRT2);

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, round_shaper);
  vect3d bf;
  dyn.getBField (vect3d(-rr, 0, 0), bf);
  bf.print ("bf at round coil center");
  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.25, 0.25*HEIGHT/WIDTH);
  screen.set_axis (vect3d(-radius * 1.15, radius * 1.15, 0), vect3d(1, 0, 0), vect3d(0, 1, 0));
  scene.calc (screen);
  fmin = scene.get_fmin();
  fmax = scene.get_fmax();
  scene.render_map (screen);
  vect3d p0, p1;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01);
  scene.find_untouching_bline (p1, c2, mid, 1, 0.01);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_round_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_round_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_round_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_round_bfield_lic_logcolor.bmp");
  }

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, round_shaper);
  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.75, 0.75*HEIGHT/WIDTH);
  scene.calc (screen);
  scene.calibrate_field (fmin, fmax);
  scene.render_map (screen);
  vect3d p0, p1;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01);
  scene.find_untouching_bline (p1, c2, mid, 1, 0.01);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_overview_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_overview_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_overview_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_overview_bfield_lic_logcolor.bmp");
  }

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, rect_shaper);
#if 0
      FILE *fout = fopen ("polywell.inc", "w");
      fprintf (fout, "%s", statics.get_povray_script().c_str());
      fclose (fout);
      exit(0);
#endif
  vect3d bf;
  dyn.getBField (vect3d(-rr, 0, 0), bf);
  bf.print ("bf at rect coil center");

  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.25, 0.25*HEIGHT/WIDTH);
  screen.set_axis (vect3d(-radius * 1.15, radius * 1.15, 0), vect3d(1, 0, 0), vect3d(0, 1, 0));
  scene.calc (screen);
  scene.calibrate_field (fmin, fmax);
  scene.render_map (screen);

  vect3d p0, p1;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01);
  scene.find_untouching_bline (p1, c2, mid, 1, 0.01);

  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_rect_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_rect_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_rect_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  scene.render_bline (screen, p1, Screen::BLUE, 1, p1, 0.01, false);
  screen.write ("polywell_cube_rvr_rect_bfield_lic_logcolor.bmp");
  }

  c1 = vect3d(-radius / PREC_SQRT2, -radius / PREC_SQRT2, rr);
  c2 = vect3d(-rr, -radius / PREC_SQRT2, radius / PREC_SQRT2);
  c3 = vect3d(-radius / PREC_SQRT2, -rr, radius / PREC_SQRT2);
  mid = (c1 + c2 + c3) / 3.0;

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, round_shaper);

  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.25, 0.25*HEIGHT/WIDTH);
  screen.set_axis (c1, vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));

  scene.calc (screen);
  scene.calibrate_field (fmin, fmax);
  scene.render_map (screen);

  vect3d p0;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01); // we calculate these only for the big round one

  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_round_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_round_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_round_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_round_bfield_lic_logcolor.bmp");
  }

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, round_shaper);

  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.8, 0.8*HEIGHT/WIDTH);
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));

  scene.calc (screen);
  scene.calibrate_field (fmin, fmax);
  scene.render_map (screen);

  vect3d p0;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01); // we calculate these only for the big round one

  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_overview_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_overview_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_overview_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_overview_bfield_lic_logcolor.bmp");
  }

  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  make_polywell_cube (statics, radius, wr, spacing, I, 0, rect_shaper);
  Scene scene(dyn);
  Screen screen(WIDTH, HEIGHT, 0.25, 0.25*HEIGHT/WIDTH);
  screen.set_axis (c1, vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen);
  scene.calibrate_field (fmin, fmax);
  scene.render_map (screen);
  vect3d p0;
  scene.find_untouching_bline (p0, c1, mid, 1, 0.01);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_rect_bfield.bmp");
  scene.set_coloring (32, true);
  scene.render_map (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_rect_bfield_log.bmp");
  scene.set_coloring (1);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_rect_bfield_lic.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  scene.render_bline (screen, p0, Screen::BLUE, 1, p0, 0.01, false);
  screen.write ("polywell_cube_rvr_pcp_rect_bfield_lic_logcolor.bmp");
  }

  return 0;
}

